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I.  EXECUTIVE  SUMMARY 


!» 


The  VSC  research  program  at  Sierra  Geophysics  during  FY81  was 
sub-divided  into  three  distinct  tasks.  These  tasks,  summarized  below, 
have  the  common  link  of  extending  our  understanding  of  body  wave 
propagation  from  underground  explosions  and  thus  enhancing  the  utility 
of  those  seismic  phases  in  the  yield  estimation  procedure. 

The  first  task  of  this  past  year's  program  was  the  determination  of 
receiver  functions  at  a  number  of  seismic  stations  for  back-azimuths 
appropriate  to  the  Nevada  Test  Site.  While  most  of  the  seismic  stations 
of  interest  are  not  well  situated  for  the  purpose  of  monitoring  NTS,  an 
attempt  to  extend  the  receiver  function  methodology  to  NTS  azimuths 
represents  an  important  test  of  the  analysis  technique.  This  research 
task  was  presented  in  detail  in  the  Sierra  Geophysics  Semi-Annual 
Technical  Report  for  the  period  December  1,  1980  to  May  31,  1981 
submitted  to  AFTAC/VSC  [VSC-TR-81-26  (SGI-R-81-041)] .  A  summary 
of  that  research  is  contained  in  Section  III  of  this  report. 

The  second  task  of  the  FY81  program  was  an  attempt  to  model  the 
effects  of  near-source  structural  effects  on  body  waves  from  under¬ 
ground  explosions.  These  effects  at  NTS  had  been  examined  previously 
and  modeling  attempted  (Hart  et  a[. ,  1979)  using  preliminary  structural 
information  and  an  optics  approach  to  propagation.  In  this  program, 
more  reliable  structural  models  were  available  as  well  as  a  more 
sophisticated  Kirchhoff  integral  equation  method.  Two  source  regions 
of  interest  were  examined.  The  results  of  modeling  propagation  in  the 
first  area,  Yucca  Flats  at  NTS,  are  presented  in  detail  in  Section  II  of 
this  report.  Also  presented  is  a  summary  of  the  theoretical  background 
of  the  Kirchhoff  integral  method.  The  results  of  the  application  of  this 

SGI-R-82-058 


2 


technique  to  a  second  test  area  are  presented  in  Part  li  of  this  report, 
submitted  separately. 

Finally,  the  third  phase  of  the  FY81  research  program  was  the 
initiation  of  a  long  term  project  to  develop  a  moment  tensor  analysis 
package  (called  MOTIVE)  to  jointly  invert  surface  waves  and  body 
waves.  This  package  will  be  able  to  incorporate  the  most  advanced 
path  corrections  for  both  data  types  and  is  intended  for  use  in  both 
seismic  discrimination  and  yield  estimation  procedures.  This  task  will 
also  determine  preliminary  moment  vs.  yield  relationships.  This  phase 
began  toward  the  end  of  the  third  quarter  of  FY81  and  is  still  basically  in 
its  preliminary  stages.  A  summary  of  the  approach  to  the  software 
development  is  contained  in  Section  IV  of  this  report.  A  detailed 
presentation  of  the  functional  specifications  of  the  package  may  be 
found  in  a  separate  technical  report  submitted  to  AFTAC/VSC 
[VSC-TR-81-022  (SGI-R-81-039)] . 
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II.  INVESTIGATIONS  OF  NEAR  SOURCE  EFFECTS  ON  BODY  WAVES, 

PART  I  -  YUCCA  FLATS 

1.  Introduction 

It  has  been  appreciated  for  some  time  that  observed  teleseismic 
magnitudes  of  underground  nuclear  explosions  at  Yucca  Flats  differ 
systematically  from  magnitudes  predicted  on  the  basis  of  known  yields 
(Alewine  et  aL ,  1977).  Figure  1  shows  a  contour  map  of  the  observed 
magnitude  residuals.  In  general,  contours  trend  north-south,  parallel 
to  the  strike  of  the  Yucca  basin,  with  anomalously  large  observed  m^ 
occurring  for  events  in  the  western  portion  of  the  basin.  This 
strongly  suggests  that  the  observed  magnitude  anomalies  are  associated 
with  shallow  basin  structure. 

There  are  two  simple  mechanisms  which  might  be  invoked  to 
explain  these  observed  anomalies.  The  first  mechanism  is  simply  a 
systematic  east-west  change  in  material  properties  at  the  working  points 
of  the  underground  explosions  in  question.  This  would  result  in 
changes  of  coupling,  and  therefore  changes  in  apparent  yield  as 
locations  change  in  the  basin.  One  might  expect,  however,  that  this 
systematic  variation  would  have  been  detected  in  the  borehole  logs. 
The  second  mechanism  is  a  focusing  and  defocusing  of  reflections  and 
transmissions  by  interfaces  which  vary  in  depth  throughout  the  basin. 
In  this  study,  we  will  restrict  ourselves  to  this  second  mechanism,  and 
examine  the  effects  of  known  basin  structures  on  both  waveform  and 
amplitude  for  a  number  of  source  sites  throughout  the  Yucca  Flats 
basin. 
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Figure  1.  (after  Alewine,  et  al.,  1977)  Residuals  of  observed  vs. 
predicted  Mw  for  sources  within  Yucca  Valley. 
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In  the  next  section  of  this  report  we  describe  the  Kirchhoff  ap¬ 
proximation,  the  specific  Boundary  Integral  Equation  method  used  in 
this  work  to  solve  the  wave  propagation  problem.  The  emphasis  of  this 
work  is  on  applications  of  this  method.  A  thorough  mathematical 
treatment  of  several  Boundary  Integral  Equation  methods,  including  the 
Kirchhoff  approximation,  is  given  elsewhere  (Wyss  and  Mellman,  1982), 
and  thus  an  heuristic,  rather  than  mathematical,  approach  is  used  in 
this  discussion.  The  aim  of  this  section  is  to  provide  the  reader  with  a 
simple  description  of  the  methods  used  in  this  work  and  the 
applicability  of  these  methods,  rather  than  providing  a  detailed 
mathematical  treatment. 

In  the  third- section  we  examine  the  synthetic  seismograms  pro¬ 
duced  using  the  Kirchhoff  approximation  on  existing  Yucca  Basin 
models,  and  compare  these  synthetic  seismograms  to  observed  data.  No 
attempt  has  been  made  to  model  details  of  specific  events,  as  this  would 
have  involved  detailed  source  modeling  far  beyond  the  scope  of  the 
current  project.  In  addition,  uncertainties  in  available  basin  structures 
and  uncertainties  in  the  effects  of  near-source  non-linearities  make  the 
results  of  such  detailed  studies  uncertain,  at  best.  Instead,  we 
examine  only  the  broad  trends  exhibited  by  the  data.  These  trends 
take  the  form  of  a  systematic  east  to  west  magnitude  bias  and  an  east 
to  west  evolution  of  wave  shape.  A  comparison  of  synthetic  seismo¬ 
grams  with  data  indicate  that  known  structure  is  sufficient  to  explain 
both  magnitude  and  waveform  changes  in  the  broad  sense. 
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2.  The  Kirchhoff  Approximation 

The  efficient  numerical  propagation  of  waves  through  irregular 
layers  of  constant  material  in  two  and  three  dimensions  has  been  a 
problem  of  considerable  geophysical  interest  of  the  last  few  years,  and 
has  received  a  great  deal  of  theoretical  consideration.  However,  even 
for  the  case  of  acoustic  wave  propagation,  this  problem  has  been  par¬ 
ticularly  difficult  numerically.  This  is  primarily  due  to  the  fact  that 
although  the  equations  of  motion  are  linear  in  the  field  quantities  of 
displacement  and  velocity,  they  are  non-linear  in  terms  of  interface 
structure.  This  fundamental  non-linearity  precludes  the  construction  of 
solutions  for  complex  structures  by  superposition  of  the  solutions  for 
simple  structures,  and  forces  one  into  computationally  costly  schemes. 

Techniques  for  dealing  with  this  difficulty  have  run  the  gamut 
from  classical  ray  techniques  to  finite  difference  codes,  with  no  single 
technique  proving  entirely  satisfactory.  In  cases  where  significant 
diffraction  and  interference  effects  are  present,  an  "exact"  solution  is 
required  and  finite  difference  methods  have  generally  been  employed. 
Unfortunately,  these  computational  techniques  are  well  known  for  their 
heavy  computational  demands  even  in  two  dimensions,  and  their  seem¬ 
ingly  inevitable  use  in  certain  problems  has  forced  the  numerical  com¬ 
munity  into  requiring  ever  larger  and  faster  computers. 

The  large  computational  requirements  of  finite  difference  methods 
are  created  by  the  necessity  of  imposing  a  grid  on  a  large  volume  of 
space.  For  problems  involving  wave  propagation  in  irregular  layers  of 
constant  velocity  material,  however,  a  more  concise  and  efficient 
formulation  is  possible.  This  formulation  leads  to  Boundary  Integral 
Equation  (BIE)  techniques,  which  require  the  gridding  of  surfaces  of 
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velocity  discontinuity  only.  This  results  in  a  much  smaller  number  of 
total  grid  points  in  BIE  calculations,  as  opposed  to  the  corresponding 
finite  difference  calculations.  The  penalty  that  must  be  paid  for  this 
reduction  in  the  number  of  grid  points  is  that  the  matrices  which 
describe  the  equations  of  motion  and  boundary  conditions,  while  sparse 
in  finite  difference  calculations,  are  dense  for  BIE  methods.  These 
dense  matrices  are  a  significant  problem  in  the  implementation  of  BIE 
methods,  and  force  one  into  a  careful  consideration  of  how  to  deal  with 
them. 

A  useful  approximation  at  this  point  is  the  separation  of  the  BIE 
problem  into  two  parts.  The  first  part  is  the  propagation  of  the  wave 
field  from  one  discontinuity  to  an  adjacent  discontinuity.  The  second 
part  consists  of  the  interaction  of  the  wave  field  with  the  discontinuity. 

The  first,  or  propagation,  stage  may  be  considered  simply  in  terms 
of  the  propagation  problem  of  a  point  source  in  a  whole  space.  Each 
grid  point  on  the  first  boundary  represents  a  point  source  in  the  whole 
space,  and  each  grid  point  on  the  second,  or  receiving,  boundary  acts 
as  a  receiver  in  a  whole  space.  Thus,  the  response  that  the  receiver 
boundary,  ignoring  boundary  interaction,  may  be  given  exactly  by 

u(x)  =  f  (u(x  )  |g  (x, x  )  “  |S  )  G(x,x  )\dS  (1) 
J\o  o  o  o/l 

x  e  S 
o  1 

where  G  is  the  whole  space  Green's  functions  and  n  is  the  unit  normal 
to  the  surface  at  xQ.  If  the  separation  between  source  and  receiver 
surfaces  is  never  less  than  approximately  one  grid  point  spacing,  a 
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stable  and  accurate  time  domain  estimate  of  u  (x)  and  du(x)/dn  may  be 
obtained  in  a  simple  straightforward  manner. 

The  second,  or  interaction  stage,  may  be  approached  in  several 
ways.  In  this  study,  we  use  the  simplest  treatment  of  the  interaction, 
which  is  known  as  the  Kirchhoff  approximation.  In  this  approximation, 
only  the  local  values  of  the  wave  field  and  the  local  direction  of  the 
surface  normal  are  considered.  Thus,  at  each  point,  the  interaction  of 
the  incident  wave  with  the  surface  is  calculated  as  if  the  local  portion 
of  the  incident  wave  were  part  of  a  plane  wave  and  the  local  boundary 
were  in  fact  a  plane  boundary  with  normal  direction  equal  to  the  local 
boundary  normal.  The  direction  and  amplitude  of  the  equivalent 
incident  plane  wave  at  each  said  point  at  a  single  point  in  time  must  be 
calculated  from  the  values  of  du/dn  and  du/dt  of  the  incident  wave  field 
at  the  appropriate  space  and  time  point.  Boundary  values  of  the  field 
variables  may  then  be  computed  through  simple  local  use  of  plane  wave 
reflection  and  transmission  coefficients.  In  general,  these  coefficients 
will  change  as  a  function  of  both  time  and  position  on  the  surface. 

It  is  a  relatively  simple  matter  to  cascade  the  two  boundary 
procedure  and  in  this  way  obtain  propagation  through  a  stack  of 
layers.  Thus,  if  we  number  our  interfaces  1  through  N,  with  an  input 
wave  incident  on  layer  1,  we  apply  the  two  boundary  procedure  first  to 
interfaces  1  and  2,  and  then  use  the  newly  computed  boundary  values 
at  interface  2  to  propagate  the  wavefield  to  layer  3.  This  procedure  is 
continued  until  interface  N  is  reached.  In  each  case,  the  Green's 
function  used  is  the  whole  space  Green's  function  appropriate  for  the 
material  through  which  propagation  is  currently  taking  place. 
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First  order  reflection  may  be  included  in  the  solution  in  a 
relatively  simple  manner.  Beginning  at  interface  N,  we  propagate  the 
boundary  values  previously  calculated  at  interface  N  to  interface  N-1. 
The  sum  of  these  boundary  values  and  the  boundary  values  computed 
as  a  result  of  propagation  from  interface  N-2  now  form  the  new  estimate 
of  the  boundary  values  at  interface  N-1.  This  new  estimate  is  then 
used  to  propagate  the  reflected  wavefield  to  layer  N-2.  This  procedure 
is  repeated  until  interface  1  is  encountered.  At  interface  M,  we  have 
estimates  of  the  boundary  values  resulting  from  propagation  from  inter¬ 
face  M-1  and  estimates  of  the  boundary  values  resulting  from  propaga¬ 
tion  from  interface  M+1.  The  total  boundary  values  at  interface  M  are  a 
sum  of  these  quantities.  By  repeating  the  propagation  procedure  "up" 
and  "down"  the  stack  of  layers  and  updating  the  appropriate  "upgoing" 
and  "downgoing"  boundary  values  at  each  interface,  it  is  possible  to 
include  multiple  reflections  to  any  desired  order.  This  procedure  may 
be  viewed  as  a  kind  of  ray  expansion  for  BIE  methods.  It  has  the 
advantage  that  no  more  than  two  interfaces  need  be  considered  at  a 
time,  and  that  the  amount  of  computation  required  is  linearly 
proportional  to  both  the  order  number  of  the  multiple  reflections  desired 
and  the  number  of  layers  in  the  stack. 

The  Kirchhoff  approximation  bears  a  strong  resemblance,  in  many 
ways,  to  geometric  optics.  In  fact,  it  is  possible  to  derive  geometric 
optics  as  a  high  frequency  approximation  to  the  Kirchhoff  BIE  formula¬ 
tion.  Thus,  the  Kirchhoff  approximation  may  be  expected  to  be  very 
accurate,  subject  to  restrictions  imposed  by  the  finite  grid,  in  domains 
where  optics  is  valid.  Unlike  optics,  the  Kirchhoff  approximation  does 
not  show  abrupt,  discontinuous  behavior  at  shadow  zone  boundaries, 
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instead  showing  the  continuous,  frequency  dependent  behavior  one 
expects  from  diffraction  theory. 

in  general,  the  Kirchhoff  approximation  gives  very  accurate 
results  in  situations  dominated  by  kinematic  effects.  Where  dynamic 
effects  are  important,  the  results  become  increasingly  inaccurate.  By 
dynamic  effects,  we  mean  those  effects  caused  by  multiple  interactions 
of  the  wavefield  with  a  boundary.  Examples  of  such  effects  include 
head  waves  and  diffraction  of  post-critical  waves.  It  is  somewhat 
difficult  to  assess  exactly  where  the  Kirchhoff  approximation  becomes 
inaccurate,  since  arrivals  usually  have  the  correct  arrival  time, 
although  possibly  an  incorrect  amplitude.  Thus,  the  Kirchhoff  approx¬ 
imation  usually  gives  plausible,  if  not  necessarily  correct,  results.  The 
question  of  accuracy  of  various  BIE  techniques  is  considered  further  in 
Wyss  and  Mellman  (1982). 

Two  models  of  the  Yucca  Basin  have  been  considered  in  this  work. 
The  northern  model  is  smooth  and  does  not  have  radically  dipping 
structures.  It  is  felt  that  the  Kirchhoff  approximation  provides  very 
accurate  results  for  this  case.  The  southern  model  contains  faults  and 
a  steeply  dipping  sidewall.  Reflected  energy  for  sources  near  the 
sidewall  is  well  beyond  critical  angle.  For  this  reason,  the  results  of 
the  southern  model  calculations  should  be  treated  as  somewhat  suspect. 
Additional  work  with  more  accurate  BIE  methods  is  necessary  for  this 
southern  model  in  order  to  assess  the  reliability  of  our  results. 

Several  additional  approximations  and  transformations  proved 
desirable  in  the  current  study.  Since  the  problem  of  interest  involves 
numerous  sources  observed  at  relatively  few  teleseismic  receivers,  a 
reciprocal  calculation  is  much  more  efficient  than  the  solution  of  the 


true  problem.  If  G..(x,x')  is  the  displacement  in  the  i  direction  at 
point  x  caused  by  a  point  force  in  the  j  direction  at  point  x*  then,  by 
reciprocity  (Aki  and  Richards,  1979) 

Gij(x,x')  =  G..(x',x)  (2) 

Differentiating  we  find 

Gij,j(X'X’)  =  G j i , j ( * ' » x ) •  (3) 

Thus,  the  equivalent  problem  to  determining  the  vertical  displace¬ 
ment  at  x  caused  by  an  explosion  at  x1  is  to  find  the  divergence  of  the 
displacement  at  x*  caused  by  a  vertical  point  force  at  x.  Letting  the 
displacement  at  x1  be  u,  we  have,  in  terms  to  the  usual  potential 
decomposition 

u  =  +  VxA,  or 

2  1  8*4 

'■!  =  '•=?  ^2-  (4) 

Thus  we  may  find  the  vertical  displacement  at  a  distant  receiver  for  an 
explosion  source  in  a  basin  by  taking  the  second  time  derivative  of  the 
dilitational  potential  at  the  source  point  induced  by  a  vertical  point 
force  at  the  receiver  point.  For  receivers  in  the  30°-90°  distance 
range,  this  reciprocal  problem  amounts  to  determining  the  dilitational 
potential  at  all  possible  source  points  for  a  plane  wave  incident  on  the 
basin  structure.  This  is  not  the  same  as  simply  interchanging  the 
source  and  receiver  points. 
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As  only  the  dilitational  potential  contributes  to  the  final  result, 
and  since  we  will  consider  only  a  short  time  window  following  the 
P-wave  first  arrival,  the  only  shear  wave  energy  that  contributes  to 
the  final  solution  must  undergo  at  least  two  mode  conversions.  As  this 
is  not  an  efficient  process,  and  as  such  energy  will  in  any  case  be 
somewhat  late  in  the  record,  we  do  not  include  shear  conversions. 
Correct  elastic  reflection  coefficients  are  used,  however,  so  that  energy 
lost  from  P-waves  due  to  such  conversions  is  accounted  for. 

Since  we  are  modeling  only  a  restricted  portion  of  the  seismogram, 
no  attempt  has  been  made  to  suppress  arrivals  due  to  grid  edge 
effects.  Such  arrivals  are  quite  evident  and  generally  occur  several 
seconds  after  the  first  arrival. 

All  synthetic  seismograms  presented  in  this  study  include  up  to 
second  order  reflections,  which  includes  all  "rays"  undergoing  a  total  of 
five  or  fewer  reflections.  This  proved  sufficient  in  terms  of  the  time 
window  of  interest  and  constraints  imposed  by  the  arrival  of  the  edge 
"reflection". 
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3.  Yucca  Flats 

The  basin  structure  is  dominated  by  two  discontinuities.  The 
shallower  of  these,  corresponding  to  the  location  of  the  water  table,  is 
known  via  well  logs  to  occur  at  a  more  or  less  constant  depth  of  about 
500  m.  throughout  the  basin.  The  deeper  discontinuity,  which  rep¬ 
resents  the  Cenozoic-Paleozoic  basement  contact,  has  been  extensively 
studied  by  Herrin  and  Goforth  (1981),  using  a  combination  of  data  from 
three  seismic  lines  and  gravity  data  acquired  throughout  the  basin.  A 
joint  inversion  was  performed  by  Herrin  and  his  colleagues  which,  in 
effect,  used  the  gravity  data  to  interpolate  between  the  seismic  data. 

A  contour  map  of  the  resulting  model  may  be  seen  in  Figure  2. 

i 

The  northern  portion  of  the  Herrin  model  in  Figure  2  is  suitable 
for  use,  with  minor  modification,  in  our  basin  model.  This  is  not  true 
of  the  southern  portion,  since  the  abrupt  eastern  sidewall  is  apparently 
caused  by  a  lack  of  data,  rather  than  representing  an  actual  feature  of 
the  basin.  As  the  steepness  of  the  basin  sidewalls  has  a  major  effect 
on  both  amplitudes  and  waveforms  for  sources  throughout  the  basin,  it 
was  decided  to  use  a  two-dimensional  model  based  on  Herrin's  profile 
E-3,  shown  in  Figure  3,  to  represent  the  southern  portion  of  the  basin, 
rather  than  the  inversion  result.  An  examination  of  Figure  2  shows 
that  the  two  dimensional  approximation  used  for  the  southern  portion  of 
the  basin  is  well  justified  by  the  three-dimensional  model. 

We  will  first  consider  sources  in  the  northern  portion  of  the  basin. 
The  model  used  is  essentially  the  Herrin  and  Goforth  model  of  the 
Paleozoic-Cenozoic  contact  with  slight  changes  necessitated  by  the 
numerical  requirements  of  the  Kirchhoff  integral  code.  Numerical  in¬ 
stabilities  due  to  inadequate  surface  sampling  can  result  if  the  distance 
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Figure  2.  Contour  nap  of  Cenozoic-Paleozoic  contact  for  Yucca 
Flats  as  determined  by  Herrin  et  al.  Line  A-A' 
corresponds  to  cross  section  shown  in  modeling  plots. 
Locations  1-9  are  selected  events  used  in  waveform 
comparisons . 
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between  boundaries  becomes  appreciably  less  than  the  surface  sampling 
distance.  For  this  reason,  the  sides  of  the  basin  were  tapered  at  a 
depth  of  200  m.,  rather  than  allowing  the  tuff  layer  to  pinch  out,  as 
occurs  in  the  Herrin  model.  Also,  in  order  to  minimize  the  effects  of 
the  boundary  modification,  it  was  necessary  to  exclude  the  water  table 
interface  from  the  model  and  use  an  average  of  saturated  and 
unsaturated  material  properties,  shown  in  Table  1,  for  the  tuff  layer. 
The  alternative  would  have  been  to  taper  the  basin  at  400  m.,  which 
would  have  had  a  more  noticeable  effect  on  the  results  than  the  course 
chosen. 

Source  locations  were  chosen  on  a  line  A-A1,  shown  on  Figure  4, 
and  were  uniformly  placed  at  a  depth  of  550  m.  The  line  A-A'  was 
chosen  to  be  perpendicular  to  the  major  axis  of  the  basin,  in  a  region 
which  is  essentially  two-dimensional  in  nature.  This  allows  for  results 
obtained  on  line  A-A*  to  be  applied  to  nearby  source  locations  with  some 
degree  of  confidence.  A  cross-section  view,  along  line  A-A1,  of  the 
model  is  shown  in  Figure  5.  Squares  represent  grid  points,  while  the 
normals  to  the  surface  at  each  grid  point  designated  by  a  short  line 
segment  at  each  grid  point.  Source  locations  are  shown  as  an  inter¬ 
mediate  layer,  and  are  numbered  from  east  to  west.  We  note  that  even 
grid  point  spacing  is  not  a  requirement  of  the  program,  and  that 
irregular  spacing  has  been  employed  to  provide  improved  accuracy  and 
efficiency. 

Since  seismograms  are  to  be  computea  for  a  number  of  source 
locations  but  only  relatively  few  receiver  locations,  computational  effort 
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TABLE  1 


Material  Properties  for  Yucca  Flats  Models 


Northern  Model 


Layer  No. 

VP 

Vs 

£ 

1 

2.4 

1.4 

2.2 

2 

4.8 

2.8 

2.6 

Southern  Model 


Layer  No. 

VP 

Vs 

£ 

1 

2.1 

1.2 

1.8 

2 

3.0 

1.8 

2.2 

3 

4.8 

2.8 

2.6 
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Figure  4.  Expanded  version  of  Figure  2  showing  contour  map  of  the 
Cenozoic-Paleozoic  contact.  Numbers  on  line  A-A' 
indicate  source  locations  for  synthetic  seismograms. 
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can  be  greatly  reduced  by  employing  a  reciprocal  calculation  scheme. 
Incident  plane  waves,  with  ray  parameter  p  =  .07  sec/km, 
corresponding  to  a  distance  of  about  60  degrees  were  used  for  the  in¬ 
coming  energy.  Three  azimuths  were  utilized,  and  have  for  simplicity 
been  labeled  west,  north,  and  east,  although  the  true  directions  are 
parallel  and  perpendicular  to  the  line  A-A1  and  thus  are  oriented  about 
20  degrees  from  the  nominal  compass  directions.  In  all  cases,  reflect¬ 
ions  up  to  order  two,  which  may  include  as  many  as  five  reflections, 
have  been  included. 

For  each  azimuth,  we  present  a  set  of  four  figures.  The  first  of 
these  shows  the  model,  with  an  arrow  indicating  the  direction  of  prop¬ 
agation  of  the  reciprocal  plane  wave.  The  second  figure  shows  the  step 
function  response  at  the  far  field  station  for  every  fifth  source  location 
on  line  A-A'.  The  third  figure  is  a  synthetic  seismogram  produced  by 
convolving  the  time  derivative  of  the  step  function  response  shown  in 
the  second  figure  with  a  Ricker  wavelet  with  time  constant  of  one 
second.  The  Ricker  wavelet  has  the  form  d/dt(e  ), 
and  provides  a  simple  means  of  producing  realistic  synthetic 
seismograms.  The  fourth  figure  in  each  set  contains  true  synthetic 
seismograms,  consisting  of  the  time  derivative  of  the  step  function 
responses  (shown  in  the  second  figure)  convolved  with  a  von  Seggern 
and  Blandford  (1972)  time  function  with  k  =  10  and  B  =  2,  an 
attenuation  operator  (Futterman,  1966)  characterized  by  t*  =  .6,  a 
typical  relatively  simple  receiver  function  and  a  KS  36000  instrument 
response. 
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The  western  azimuth  is  shown  in  Figures  5  through  8.  As 
expected,  the  step  responses  in  Figure  6  show  delays  in  arrival  time  as 
the  source  moves  to  the  east.  The  arrival  from  source  location  1 
appears  somewhat  larger  than  the  arrival  from  source  6,  probably  due 
to  focusing  of  the  direct  ray  by  the  western  portion  basin.  There  is 
also  a  small  change  in  pP  time  between  sources  1  and  6,  which  is 
denoted  by  a  return  to  zero  of  the  step  response  due,  as  before,  to 
the  western  sidewall.  The  most  significant  feature,  however,  is  the 
change  in  arrival  time  and  amplitude  of  the  second  negative  phase, 
which  represents  a  large  reflection  off  the  bottom  and  side  of  the 
basin.  Due  to  focusing,  this  arrival  is  very  large  and  early  for  source 
1,  becoming  progressively  smaller  and  more  delayed  for  sources  toward 
the  center  of  the  basin  (e.g.,  source  16).  This  arrival  then  increases 
in  amplitude  and  its  delay  lessens  for  sources  near  the  eastern  side  of 
the  basin,  a  result  of  focusing  of  the  reflection  off  the  eastern  wall. 

A  second  order  positive  reflection  may  also  be  seen,  particularly 
for  source  16,  although  this  phase  arrives  too  late  to  have  an  effect  on 
the  early  portion  of  the  waveform  (and,  hence,  on  mfa).  The  final 
arrival,  a  ramp  function,  is  a  truncation  phase  caused  by  "reflections" 
off  the  edge  of  the  model.  Fortunately,  this  phase  arrives  with 
sufficiently  large  time  delay  and  with  sufficiently  long  period  that  it 
does  not  significantly  affect  the  synthetic  seismograms. 

The  influences  of  the  multiple  reflections  are  quite  obvious  in 
Figures  7  and  8.  The  large,  early  reflection  for  source  1  interferes 
constructively  and  produces  a  very  large,  narrow  first  trough  and 
second  peak  in  the  waveform.  The  greater  delay  for  source  6  results 
in  broadening  of  the  waveform  in  general  and  an  attendant  decrease  in 
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Figure  6.  Step  function  response  for  teleseismic  receiver  (p 
.07  sec/km)  west  of  basin.  Number  of  response 
corresponds  to  source  location  number  in  Figure  5. 
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Figure  8.  Synthetic  seismogram*  for  model  in  Figure  5  using 
SRO  instrument,  von  Seggern-Blandford  source 
function  with  k  ■  10,  B  ■  2,  receiver  function 
for  station  E  from  E.  Kazakh  events  and  additional 
attenuation  operator  with  t*  ■  .6  sec.  Western 
azimuth. 
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Figure  10.  Step  function  responses  for  North  azimuth 
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Figure  14.  Step  function  responses  for  East  azimuth 
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amplitude.  As  the  multiple  reflections  move  still  further  back  in  time, 
for  sources  11  and  16,  we  lose  the  constructive  interference  effect, 
with  a  resultant  decrease  in  amplitude  and  the  development  of  an  in¬ 
flection  in  the  second  peak.  Near  the  eastern  portion  of  the  basin, 
(sources  21  and  26),  we  find  an  increase  in  amplitude  and  decrease  in 
pulse  width,  due  to  the  proximity  of  the  eastern  wall. 

A  comparison  of  Figures  7  and  8  shows  that  while  details  of  wave¬ 
forms  differ,  the  overall  trends  in  amplitude  and  waveform  are  the  same 
for  the  Ricker  wavelet  and  explosion  source  synthetics.  This  fact 
indicates  that  even  though  these  effects  are  caused  by  interference  of 
several  arrivals,  they  are  not  unduly  sensitive  to  small  changes  in 
source  function. 

The  model  and  synthetic  seismograms  appropriate  to  a  northern 
azimuth  are  presented  in  Figures  9-12.  While  the  northern  azimuth 
results  differ  in  detail  from  those  for  the  eastern  azimuth  (Figures 
13-16),  the  overall  trends  are  the  same,  in  terms  of  both  amplitude  and 
waveform.  In  general,  somewhat  less  overall  amplitude  variation  is 
present  for  both  of  these  azimuths  in  comparison  to  the  amplitude 
variation  occurring  for  a  western  azimuth.  It  is  also  somewhat  sur¬ 
prising  that  synthetics  for  source  26  show  smaller  amplitudes,  at  all 
azimuths,  than  synthetics  from  other  source  locations.  This  appears  to 
be  related  to  the  large  slope  of  the  basin  to  the  east.  It  is  not  known 
at  this  time,  however,  whether  this  effect  is  real  or  merely  an  artifact 
of  the  Kirchhoff  approximation. 

The  predicted  amplitude  anomaly  as  a  function  of  source  location 
and  azimuth  is  summarized  in  Table  2.  No  attempt  has  been  made  here 
to  correct  for  either  period  or  instrument.  In  terms  of  general  trend, 
the  decrease  in  apparent  magnitude  from  the  western  edge  to  the  center 
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of  the  basin  agrees  well  with  the  results  of  Alewine  et  aL  (1977)  shown 
in  Figure  1.  Comparison  with  a  more  recent  study  by  Alewine  (per*- 
sonal  communication),  shown  in  Figure  17,  do  not,  however,  show  good 
agreement  in  detail.  For  reference,  Figure  18,  which  shows  line  A-A', 
is  scaled  identically  with  Figure  17.  It  cannot  be  presently  determined 
whether  the  discrepancy  is  due  to  lack  of  a  sufficiently  accurate  or 
detailed  basin  model,  to  inadequacies  in  m^  as  a  measure  of  event  size 
when  interference  effects  produce  large  changes  in  apparent  period,  or 
to  some  other  factor. 

The  predicted  evolution  of  waveforms  as  a  function  of  location  in 
the  basin  provides  a  test  of  the  validity  of  the  modeling  procedure. 
Figure  18  shows  the  basement  structure,  the  line  A-A  and  the  location 
of  nine  underground  nuclear  explosions  in  reasonably  close  proximity  to 
line  A-A'.  Available  seismograms  for  these  events  at  ZOBO,  MAJO  and 
NAO  are  shown  in  Figures  19-21.  A  projection  of  these  events  onto 
line  A-A1  was  done  in  such  a  way  as  to  preserve  the  distance  from  the 
eastern  edge  of  the  basin.  The  time  scale  of  the  observed  seismograms 
in  the  figures  is  the  same  scale  used  in  the  synthetics. 

Due  to  the  symmetry  of  the  basin,  seismograms  recorded  at  ZOBO 
(a  southerly  azimuth),  Figure  19,  should  be  comparable  to  synthetics 
for  a  northern  azimuth.  The  evolution  of  the  inflection  in  the  second 
peak  and  development  of  the  second  trough  as  predicted  by  the  syn¬ 
thetics  when  the  source  is  moved  toward  the  center  of  the  basin  can 
clearly  be  seen  in  the  data.  The  trend  in  the  synthetic  waveforms  is 
thus  confirmed,  despite  the  fact  that  no  attempt  was  made  to  correct 
source  time  functions  or  depths  for  individual  events. 


SGI-R-82-058 


YUCCA  VALLEY  MAGNITUDE  RESIDUALS 
(OBSERVED  -  PREDICTED) 


FARALLONES 


L / 


REBLOCHON 


MARSILLY 


X 


HEARTS 


SCANTLING 


RUDDER 


RUMMY 


SANDREEF 


0 


116°  WEST  LONGITUDE 
YUCCA  VALLEY  MAGNITUDE  RESIDUALS 
(OBSERVED  •  PREDICTED) 


Figure  17.  (after  Alewine  ,  et  al.  1977)  Observed  magnitude 
residuals  for  Yucca  Flats. 
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Figure  18.  Expanded  version  of  Figure  2  showing  depth  to 
Cenosoic-Paleozoic  contact,  in  km.  Line  A-A' 
corresponds  to  cross  sections  in  Figures  5,  9,  and 
13,  with  source  locations  indicated.  Hexagonal 
symbols  indicate  locations  of  events  used  for 
waveform  comparisons. 
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The  seismograms  for  MAJO,  shown  in  Figure  20,  show  a  great 
similarity  to  the  ZOBO  seismograms.  The  moveout  of  the  Cenozoic- 
Paleozoic  interface  reflection  is  even  more  obvious  for  this  station  than 
for  ZOBO.  Once  again,  good  qualitative  agreement  is  obtained  between 
synthetic  and  observed  waveforms,  although  it  appears  that  the  syn¬ 
thetics  have  somewhat  underestimated  the  size  of  the  reflected  arrival. 

Seismograms  recorded  at  NAO  once  again  exhibit  approximately  the 
same  predominant  period  as  do  the  synthetics.  As  with  the  other  data, 
we  see  a  general  trend  toward  broadening  and  development  of  an  in¬ 
flection  as  we  move  the  source  toward  the  center  of  the  basin.  For 
this  station,  however,  it  is  more  apparent  than  for  the  previous  two 
stations,  that  detailed  agreement  between  data  and  synthetics  is  lacking 
although  general  trends  do  appear  to  be  correct. 

We  now  consider  sources  in  the  southern  portion  of  the  basin.  A 
two-dimensional  model,  based  on  the  Herrin  E-3  profile  (shown  in 
Figure  3  with  material  properties  shown  in  Table  1)  was  used  in  this 
phase  of  the  study.  It  was  also  necessary  to  modify  this  model  to 
insure  numerical  stability  of  the  Kirchhoff  integral  code.  This  was 
done  by  tapering  both  ends  of  the  basin  in  order  that  no  layer  thick¬ 
ness  is  reduced  below  200  m.  We  have  elected,  in  this  case,  to  retain 
the  water  table  in  the  model.  The  intermediate  interface  plays  the  dual 
role  of  source  location  and  boundary  between  saturated  and  unsaturated 
tuff.  Sources  are  located  just  below  the  water  table. 
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Figure  19. 


Waveforms  of  events  proceeding  east  to  vest  across 
basin,  recorded  at  station  ZOBO.  Note  development 

of  interference  effect  on  back  of  second  peak. 
Traces  scaled  to  equal  amplitude  by  plotting 
routine . 
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Figure  21.  Waveform*  recorded  at  station  nao  for  events  proceeding 

west  to  east  across  basin. 
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Models  and  results  for  the  northern  azimuth  are  presented  in 
Figures  22-25  and  for  the  eastern  azimuth  in  Figures  26-29.  The 
synthetic  seismograms  in  Figures  24,  25,  28  and  29  show  a  large 
amplitude  increase  for  both  azimuths  for  source  location  7,  located  very 
near  the  boundary  fault.  Moving  the  source  further  east  produces  a 
relatively  slight  amplitude  decrease  with  minimum  amplitude  occurring  at 
source  location  16,  followed  by  an  amplitude  increase  as  the  source  is 
moved  still  further  east.  Some  waveform  variation  is  observed  in  the 
eastern  part  of  the  basin  but,  with  the  exception  of  source  location  7, 
very  little  waveform  change  is  evident  in  the  western  portion. 
Comparison  of  the  amplitude  behavior  of  the  synthetics  with  the 
observed  magnitude  residuals  in  Figures  1  and  3  show  some  qualitative 
agreement. 

Both  the  amplitude  and  waveform  variation  of  the  northern  basin 
model  contrast  sharply  with  that  observed  for  the  southern  basin  model. 
The  differences  appear  to  be  related  to  the  extreme  dip  of  the  basement 
near  the  west  side  of  the  basin  in  the  southern  model.  Numerical 
experiments  have  shown  that  the  dip  of  this  boundary  fault  controls  the 
amplitude  of  the  basement  reflection  which,  in  turn,  creates  the  major 
portion  of  the  amplitude  and  waveform  variation  observed  for  the 
northern  model.  As  the  dip  of  this  feature  increases,  the  amplitude  of 
the  reflection  goes  through  a  maximum  for  source  locations  east  of  the 
fault. 

Since  the  Kirchhoff  approximation  does  not  correctly  handle  critical 
reflections  and  most  diffractive  effects,  both  of  which  will  be  present 
when  extreme  dips  in  structure  occur,  these  observed  effects  could 
conceivably  be  an  artifact  of  the  modeling  technique.  Further  study 
would  be  required  to  answer  this  question. 
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Figure  22.  Herrin  et  al.  profile  E-3,  used  as  cross  section  for  two-dimensional  southern 
basin  model.  Grid  shown  for  surface,  water  table  and  Paleozoic  interface. 

Sources  located  directly  below  water  table.  For  receiver  at  northern  or 
southern  azimuth. 
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Figure  23.  Step  function  responses  for  model  in  Figure  23  at  north 


or  south  azimuth. 
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Figure  25.  Synthetic  seisroogramB  for  explosion  sources,  as  in 
Figure  8,  but  for  model  in  Figure  22  with  receiver 
at  north  or  south  azimuth. 
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Model  as  in  Figure  22  but  for  receiver  at  east  azimuth 
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Figure  27.  Step  function  responses  for  sources  in  Figure  26 
Receiver  at  east  azimuth. 
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Figure  28.  Ricker  wavelet  synthetic  seismograms  for  model  in 
Figure  26.  Receiver  is  at  east  azimuth  with 
p  *  .07  sec  Am. 


1147 


4.  Conclusions  and  Recommendations 


The  Kirchhoff  integral  approximation  has  proved  to  be  a  useful 
tool  in  examining  the  effects  of  shallow  structure  on  body  wave  amp¬ 
litudes  and  waveforms  from  underground  nuclear  explosions.  In  the 
example  of  Yucca  Flats,  it  has  been  possible  to  broadly  model  the 
observed  waveform  variation  and  essentially  all  of  the  observed  amp¬ 
litude  variations  in  the  northern  portion  of  the  basin  solely  by  modeling 
interactions  with  the  Paleozoic-Cenozoic  contact.  Differences  in  the 
finer  details  do,  however,  exist  between  the  observations  and  pre¬ 
dictions.  Several  important  conclusions  and  observations  can  be  drawn 
from  this  modeling  effort.  Among  the  more  significant  is  the  fact  that 
the  waveform  and  amplitude  variations  are  the  result  of  interference 
effects  resulting  from  known  shallow  structure  and  are  not  the  result  of 
focusing  caused  by  deep  structures.  In  addition,  it  is  clear  that  these 
interference  effects  are  controlled  by  the  impedance  contrast  at  the 
Paleozoic-Cenozoic  contact  and,  particularly  with  respect  to  the  extent 
and  spatial  location  of  the  interference,  by  the  steeply  dipping  sidewalls 
of  the  basin.  These  facts  suggest  that  an  accurate  knowledge  of  depth 
of  burial  and  source  time  function  is  essential  to  matching  the  finer 
details  of  the  seismic  radiation  from  this  region. 

Although  one  can  confidently  state  that  the  amplitude  and  wave¬ 
form  variations  at  Yucca  Flats  have  been  satisfactorily  modeled,  there 
are  other  lines  of  related  research  that  should  be  considered  for 
further  study.  The  first  of  these  involves  the  use  of  more 
sophisticated  global  approximation  boundary  integral  equation  techniques 
to  deal  with  problems  created  by  pinch-out  of  layers  and  by  steeply 
dipping  structures.  These  are  apparent,  for  example,  near  the 
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boundary  fault  in  the  southern  portions  of  Yucca  Flats.  This  research 
is  not  required  to  confirm  our  conclusions  with  regard  to  the  Yucca 
Flats  anomalies  but  would  provide  valuable  parameters  and  under¬ 
standing  for  analysis  of  future  source  locations  of  interest. 
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III.  RECEIVER  FUNCTIONS  APPROPRIATE  TO  NTS 


The  purpose  of  this  study  is  to  estimate  relative  receiver  functions 
for  the  Nevada  Test  Site  Global  Array  (NTSGA)  which  is  composed  of 
all  high  quality  digital  seismic  stations  in  the  distance  range  30°  to  90° 
from  the  source  location.  Relative  receiver  functions  (RRF's)  model  all 
near-receiver  contributions  to  observed  waveforms  at  stations  in  the 
array  relative  to  the  reference  station.  The  computed  receiver  func¬ 
tions  will  permit  more  complete  forward  modeling  of  synthetic  seismo¬ 
grams  of  NTS  nuclear  explosions  as  observed  at  stations  within  the 
array. 

The  calculation  of  relative  receiver  functions  is  done  in  two  steps. 
First,  trace  deconvolutions  are  computed  by  spectral  ratio  of  the  trace 
at  a  secondary  station  to  the  trace  at  a  reference  station.  The  trace 
deconvolutions  are  computed  for  each  event  in  the  data  set  for  which 
both  reference  and  secondary  station  seismograms  are  available,  and  the 
results  for  each  station  pair  are  averaged  in  the  log-frequency  domain. 
The  result  is  an  average  transfer  function  which  permits  computation  of 
a  synthetic  trace  from  the  reference  seismogram. 

The  second  step  involves  the  recovery  of  the  receiver  functions 
from  the  deconvolutions  with  the  additional  constraint  that  the  receiver 
functions  are  as  simple  or  "delta-like"  as  is  consistent  with  the  informa¬ 
tion  in  the  trace  deconvolutions.  This  can  be  done  by  finding  a  linear 
oeprator  which,  when  convolved  with  trace  deconvolutions,  will  maximize 
the  order  (or  minimize  the  entropy)  in  the  set  of  estimated  receiver 
functions  as  a  whole.  Although  the  entropy  minimum  found  is  not 
unnecessarily  unique,  tests  have  shown  that  the  significant  information 
in  a  set  of  receiver  functions  determined  by  well-constrained  trace 
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deconvolutions  is  stable  as  a  function  of  reference  station,  number  of 
secondary  stations,  and  choice  of  weighting  schemes. 

There  are  16  stations  which  qualify  for  the  NTSGA,  including 
SRO,  LRSM,  and  NORSAR  stations.  The  SRO  station  ZOBO  in  Bolivia 
was  chosen  as  the  reference  station  because  it  had  more  useful  traces 
than  any  of  the  other  stations  in  the  array,  and  because  tests  showed 
that  ZOBO  seismograms  have  high  frequency  character  relative  to  the 
rest  of  the  array. 

Selecting  ZOBO  as  the  reference  station  allowed  the  use  of  17  NTS 
events  between  May  1977  and  September  1980.  Eleven  of  the  15 
secondary  stations  in  the  NTSGA  had  four  or  more  usable  trace  pairs 
with  ZOBO  for  these  events.  Although  we  have  obtained  all  useful  data 
from  the  SRO  stations  for  this  time  period,  gaps  in  data  exist  for  the 
AEDS  stations  from  May  1977  to  December  1977  and  September  1978  to 
June  1979. 

Estimated  relative  receiver  functions  were  computed  for  the  twelve 
stations  using  MED  processing.  Based  on  the  complexity  and  frequency 
content  of  the  computed  receiver  functions,  it  appears  that  ZOBO  was 
not  the  best  choice  for  a  reference  station.  However,  KONO,  which 
might  have  been  considered  as  a  better  reference  station,  failed  to 
provide  enough  data  to  be  used  as  the  reference  station. 

The  use  of  receiver  functions  as  deconvolution  filters  is  briefly 
examined.  The  problems  include  bandwidth  and  noise,  but  the  most 
important  concept  developed  is  that  receiver  functions  are  designed  as 
convolution  filters  and  thus  are  not  appropriately  used  in  a 
deconvolution  sense.  Further  research  will  be  required  to  define  filters 
which  may  be  used  to  simplify  observed  seismograms  prior  to  magnitude 
estimation . 
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IV.  MOTIVE  FUNCTIONAL  SPECIFICATIONS 

A  number  of  techniques  are  currently  available  for  seismic 
discrimination  and  yield  estimation.  Most  of  these  methods,  however, 
fail  to  fully  utilize  the  available  azimuthal  information  and,  indeed,  many 
methods  for  magnitude  estimates  use  only  a  single  amplitude  measure. 
Shallow  earthquakes  will,  in  general,  produce  waveforms  and  amplitude 
patterns  that  vary  azimuthally  while  a  pure  explosion  in  theory 
produces  azimuthally  independent  waveforms  and  amplitudes.  This  fact 
is  not  utilized  in  most  discrimination  schemes,  in  the  yield  estimation 
context,  azimuthal  variations  in  waveform  and  amplitude  may  contain 
information  on  tectonic  release  for  which  a  correction  should  be,  but  is 
usually  not,  made>  in  obtaining  yield  estimates.  In  the  technical  report 
on  the  subject  (VSC-TR-?1-22),  we  describe  a  computer  program, 
MOTIVE,  which  makes  use  of  this  azimuthal  information,  as  well  as 
utilizing  significant  portions  of  the  waveform  information,  in  both 
discrimination  and  yield  estimation  contexts. 

MOTIVE  is  a  joint  body  wave  and  surface  wave  moment  tensor 
inversion  program.  Given  properly  windowed  body  wave  and  surface 
wave  seismograms  and  a  trial  depth,  MOTIVE  determines  the  second 
order  moment  tensor  that  provides  the  best  fit,  in  a  least  squares 
sense,  of  synthetics  to  the  data.  By  using  a  number  of  trial  depths 
and  selecting  the  depth  that  produces  the  minimum  error,  depth,  as 
well  as  source  type  and  orientation,  may  be  determined. 

MOTIVE  incorporates  several  unique  features  in  addition  to  the 
joint  use  of  body  wave  and  surface  wave  data.  The  user  can  specify 
any  of  four  source  types;  an  unconstrained  moment  tensor,  an  isotropic 
source,  a  double  couple  source  and  a  double  couple  plus  isotropic 
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source.  This  gives  the  user  the  ability  to  explore  the  sensitivity  of  his 
solution  to  the  most  commonly  made  physical  assumptions  about  the 
source.  Moreover,  MOTIVE  can  include  state-of-the-art  body  wave  and 
surface  wave  path  and  receiver  corrections,  determined  for  each  source- 
station  pair.  This  should  help  eliminate  one  of  the  major  sources  of 
error  and  bias  in  moment  tensor  estimation,  particularly  where  short- 
period  body  wave  seismograms  are  used. 

A  flexible  modular  approach  has  been  used  in  developing  MOTIVE, 
in  order  that  the  same  program,  through  choice  of  user  options,  may  be 
used  for  both  routine  processing  and  research.  This  flexibility  also 
insures  that  future  improvements  can  be  readily  incorporated  into  the 
existing  program  structure. 
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Seismograms  recorded  at  NAO  once  again  exhibit  approximately  the 
same  predominant  period  as  do  the  synthetics.  As  with  the  other  data, 
we  see  a  general  trend  toward  broadening  and  development  of  an  in¬ 
flection  as  we  move  the  source  toward  the  center  of  the  basin.  For 
this  station,  however,  it  is  more  apparent  than  for  the  previous  two 
stations,  that  detailed  agreement  between  data  and  synthetics  is  lacking 
although  general  trends  do  appear  to  be  correct. 

We  now  consider  sources  in  the  southern  portion  of  the  basin.  A 
two-dimensional  model,  based  on  the  Herrin  E-3  profile  (shown  in 
Figure  3  with  material  properties  shown  in  Table  1)  was  used  in  this 
phase  of  the  study.  It  was  also  necessary  to  modify  this  model  to 
insure  numerical  stability  of  the  Kirchhoff  integral  code.  This  was 
done  by  tapering  both  ends  of  the  basin  in  order  that  no  layer  thick¬ 
ness  is  reduced  below  200  m.  We  have  elected,  in  this  case,  to  retain 
the  water  table  in  the  model.  The  intermediate  interface  plays  the  dual 
role  of  source  location  and  boundary  between  saturated  and  unsaturated 
tuff.  Sources  are  located  just  below  the  water  table. 
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